# load libraries ----
# reporting
library(knitr)
# visualization
library(ggplot2)
library(ggthemes)
library(patchwork)
library(viridis)
# remote sensing & GIS
library(sf)
library(raster)
library(ggspatial)
library(rnaturalearth)
library(rnaturalearthhires)
library(rnaturalearthdata)
# data wrangling
library(plyr)
library(dplyr)
library(tidyr)
library(tibble)
library(readr)
library(skimr)
library(janitor)
library(glue)
library(magrittr)
library(lubridate)
# stats
library(corrplot)
library(Hmisc)
# set other options ----
options(scipen=999)
knitr::opts_chunk$set(
tidy = FALSE,
message = FALSE,
warning = FALSE)35 Patterns of Topography & Vegetation II
Learning Objectives
After completing this tutorial you should be able to
- understand the importance of long-term monitoring stations and baseline data
- understand how spatial and temporal patterns affect ecological processes
- describe what raster data sets are and how they encode spatial information
- read spatial raster data in R, making simple maps, and extracting information for non-spatial statistical tests.
- test whether plant growth (greenness & height) is drive more by elevation, slope, or aspect.
- consider large scale patterns of relationships between topography and vegetation
Download the directory for this project here, make sure the directory is unzipped and move it to your bi328 directory. You can open the Rproj for this module either by double clicking on it which will launch Rstudio or by opening Rstudio and then using File > Open Project or by clicking on the Rproject icon in the top right of your program window and selecting Open Project.
There should be a file named 35_vegetation-I.qmd in that project directory. Use that file to work through this tutorial - you will hand in your rendered (“knitted”) quarto file as your homework assignment. So, first thing in the YAML header, change the author to your name. You will use this quarto document to record your answers. Remember to use comments to annotate your code; at minimum you should have one comment per code set1 you may of course add as many comments as you need to be able to recall what you did. Similarly, take notes in the document as we discuss discussion/reflection questions but make sure that you go back and clean them up for “public consumption”.
1 You should do this whether you are adding code yourself or using code from our manual, even if it isn’t commented in the manual… especially when the code is already included for you, add comments to describe how the function works/what it does as we introduce it during the participatory coding session so you can refer back to it.
Before you get started, let’s make sure to read in all the packages we will need.
35.1 Comparison of relationship of vegetation and topography across the US
Let’s take a look at the different NEON sites across the US and choose locations from the contiguous 48.
We are going to come up with a sampling design to compare impact of topography on vegetation across domains and regions. We can access the elevation and NDVI data sets from the Data Portal. First, we will want to know what data is available and then based off of that information you will be able to identify six sites in addition to SOAP that you want to explore.
Here are the instructions on how to find the NEON data products that contain the data we are interested in (dtm and dsm models and ndvi data)2.
2 Reminder: We have talked a lot about the importance of matching data sets for spatial and temporal resolution. Keep this in mind as you design your data set.
- Go to the NEON Data Portal.
- In the search bar type in “DTM”, select
Elevation - LiDARwhich contains the DTM and DSM from the search results. - This will take you to the Elevation - LiDar product (NEON Product ID DP3.30024.001) page.
- Scroll through the page to view the abstract, citation, study description and other documentation to the
Availability and Downloadsection which shows you for which sites (coded asSTATE-DOMAIN-SITEID) and for which years there are data available. - When you are ready to download data, click on the
Download Databutton.- First, indicate which sites and which date range you want to download.
- Then, click on the
Provisional Databutton and selectExclude. - Next, click on
Filesbutton and select the appropriateNEON_D##_SITE_DP3_######_#######_DSMandNEON_D##_SITE_DP3_######_#######_DTMtiffiles. You need to make sure that the 6 and 7 digit numbers in the file name which indicate the coordinates of the northeast corner of the tile are the same. One way to do this is to use the Filter box under the Name column header. - Now select the
Documentationbutton and selectExclude. - Finally, on the
Policiespage check theI agreebox. - Et voila, you can click the blue
Download Databutton to get your*.tiffiles.
- Once you download your
*.tiffiles you will want to make sure to move them to yourdatafolder in your project directory. If you are downloading multiple files at once you will need to make sure that you unzip the folder and pull out only the_DTM.tifand_DSM.tiffiles that you need. - Repeat the process to view and access NDVI data by typing
ndviin the search bar and selecting theVegetation indices - spectrometer - mosaicresult. - This will take you to the
Vegetation indices - spectrometer - mosaicdata product page (NEON Product ID is DP3.00026.001) where you can access theAvailability and Downloadsection to view and access available data - When you are ready to download data, click on the
Download Databutton.- First, indicate which sites and which date range you want to download.
- Then, click on the
Provisional Databutton and selectExclude. - Next, click on
Filesbutton and select the appropriateNEON_D##_SITE_DP3_######_#######_VegIndices.tiffiles - you must make sure to match the 6 and 7 digit numbers (coordinates of the tile) to your DTM and DSM files! - one way to do this is to type that part of the file name into the filter box of the name column header3. - Now select the
Documentationbutton and selectExclude. - Finally, on the
Policiespage check theI agreebox. - Et voila, you can click the blue
Download Databutton to get your*.tiffiles.
- Unfortunately, you will likely end up with a nested zip folder - make sure you diligently unzip until you find your
*.tiffile with the appropriate ndvi data.
3 This is where having thought about matching spatial and temporal resolution comes into play - you need to have data from the same year and the same location (coordinates).
Make sure you have matching sets of *.tif files for your DTM, DSM, and ndvi data. This means they should be from the same year and location. You can match location using the Northing/Easting coordinates that are part of the filename. When you are downloading files in the portal use the Filter box undernearth the name column name. You frequently have more hits than you can view on one screen so just scrolling through will make it difficult to find the match.
We will compile all of the results into a map, so you can see how the impact of topographic variables varies across our studies locations.
First, let’s read in a text file that has all the NEON sites and information on their geographic location. We are only interested in sites from the contiguous 48 so we will remove locations in Alaska, Hawaii, and Puerto Rico.
sites <- read_delim("data/neon_sites.txt", delim = "\t") %>%
clean_names() %>%
separate(site_type, into = c("tmp", "site_cat"), sep = " ", remove = FALSE) %>%
dplyr::select(-tmp) %>%
filter(!state %in% c("AK", "HI", "PR"))
head(sites)# A tibble: 6 × 12
site_name site_id domain_number state latitude longitude site_type site_cat
<chr> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
1 Abby Road ABBY D16 WA 45.8 -122. Relocata… Terrest…
2 Arikaree Ri… ARIK D10 CO 39.8 -102. Core Aqu… Aquatic
3 Ordway-Swis… BARC D03 FL 29.7 -82.0 Core Aqu… Aquatic
4 Bartlett Ex… BART D01 NH 44.1 -71.3 Relocata… Terrest…
5 Upper Big C… BIGC D17 CA 37.1 -119. Relocata… Aquatic
6 Blandy Expe… BLAN D02 VA 39.1 -78.1 Relocata… Terrest…
# ℹ 4 more variables: site_subtype <chr>, utm_zone <dbl>, easting <dbl>,
# northing <dbl>
Let’s start simple by creating a map that has all of our NEON sites on it.
First, we need to get a shapefile that we will use as our basemap - because we are going to want to orient ourselves within the US, we will also pull information we can use to plot state lines.
world <- ne_countries(scale = "medium", returnclass = "sf")
us_states <- ne_states(country = "United States of America", returnclass = "sf")We are not going to want to plot the entire world, so let’s identify the lat/longs for the four corners of the box we want to plot but pulling information on the minimum and maximum latitude & longitude of our NEON sites and then adding a small buffer around it.
Then we can plot our map.
# lat/long for map extend
x_min <- min(sites$longitude) - 2
x_max <- max(sites$longitude) + 2
y_min <- min(sites$latitude) - 2
y_max <- max(sites$latitude) + 2
# set font sizes for text labels
site_size <- 2
# set color for fill
map_color <- "#a8d481"
# create plot
ggplot() +
geom_sf(data = world, color = "black", fill = map_color) + # plot outline of countries
geom_sf(data = us_states, color = "black", fill = NA) + # plot outline of states
coord_sf(xlim = c(x_min, x_max),
ylim = c(y_min, y_max)) + # plot boundaries for map
geom_text(data = sites, aes(x = longitude, y = latitude, # add sites as labels
label = site_id, color = site_cat),
size = site_size) +
scale_color_manual(values = c("darkblue", "darkred"))NEON sites are intentionally designed so aquatic and terrestrial sites are ideally right next to each other, so those are plotting right on top of each other.
Instead of the site id’s we could also make the same map with points for locations.
# create plot
ggplot() +
geom_sf(data = world, color = "black", fill = map_color) + # plot outline of countries
geom_sf(data = us_states, color = "black", fill = NA) + # plot outline of states
coord_sf(xlim = c(x_min, x_max),
ylim = c(y_min, y_max)) + # plot boundaries for map
geom_point(data = sites, aes(x = longitude, y = latitude, # add sites as labels
fill = site_cat),
shape = 21, size = 2, alpha = 0.5) +
scale_fill_manual(values = c("darkblue", "darkred"))If we want to compare differences in the relationships of topographic variables across all the sites we plotted, we could use the same method to plot the correlation coefficients, in addition to the site labels; all we need to do is add the correlation coefficients into our sites data.frame4.
4 Once we have information for multiple sites, you will need to either read in all the text files with results and use bind_rows to create a single data.frame, or we will also allow for some outside R action, where you just cut and paste them all into a single text file and read that in.
Here is how you can read in all your files and combine them into a single data.frame. We will create an empty list and then fill it with individual data.frames (one for each site) by looping over a vectory with all the filenames. Because we names our correlations files in a very standardized way, we can specify a pattern, and R will create a vector with all those filenames in the results folder for us.
# path to results
path <- "results"
# pattern for filenames
pattern <- "correlation"
# list with all txt files in data directory
list.filenames <- list.files(path = path,
pattern = pattern)
# empty list to load files into
list.data <- list()
# loop to read in data files
for (i in list.filenames){
file <- as.character(glue("{path}/{i}"))
list.data[[i]] <- read_delim(file, delim = "\t")
}
correl_all <- ldply(list.data, data.frame)
# # write to file
# write_delim(correl_all, "results/Comp_pearson.txt", delim = "\t")For example, we could create a map for slope and NDVI like this:
# add correlation coefficients
sites_pears <- sites %>%
left_join(correl_all, by = c("site_id" = "Site")) %>%
filter(!is.na(pearson)) %>% # retain only sites with pearson coefficient estimated
filter(Param1 == "ndvi" & Param2 == "slope") # choose veg & topog relationship to map
# set font sizes for text labels
site_size <- 4
cor_size <- 5
# create plot
ggplot() +
geom_sf(data = world, color = "black", fill = map_color) + # plot outline of countries
geom_sf(data = us_states, color = "black", fill = NA) + # plot outline of states
coord_sf(xlim = c(x_min, x_max),
ylim = c(y_min, y_max)) + # plot boundaries for map
geom_text(data = sites_pears, aes(x = longitude, y = latitude, # add sites as labels
label = site_id),
size = site_size, color = "black") +
geom_text(data = sites_pears, aes(x = longitude, y = latitude, # add pearson coefficient as labels
label = round(pearson, digits = 2)),
size = cor_size, color = "darkred", nudge_y = -1.5) + # shift label below actual coordinates
theme(legend.position = "none")35.2 Acknowledgments
These activities are based on the EDDIE Remote Sensing module.5
5 Dahlin, K. M. (2021). Remote Sensing of Plants and Topography in R (Project EDDIE).